# There appears to be a problem with cum mode not including previous
# lines, and I'd rather write a python program than debug vpfit.  so
# here it is.

# Given a forest region, breaks it up into overlapping chunks, and
# fits all of them.  While fitting each chunk, includes all the lines
# from previous fit (only) as fixed.

import astro
import os
import sys

lya = 1215.67
lyb = 1025.72
zqso = 3.08
specfile = '../HE0940m1050un.fits'
inputfile = 'f26in.temp'
vsig = '3.0'
chi_thresh = '2.0'
vel=800.      # 1200 kms ~17 Ang at 5000 Ang
vel2=600.     # 200 km/s overlap between regions
SPDLGT=300000.
vc = vel/SPDLGT
vc2 = vel2/SPDLGT
wmax = lya * (1. + zqso)
wmin = lyb * (1. + zqso)
wrange = wmax - wmin
run_dir = 'run_c'
logfile = 'log.txt'

# New run, make a directory.
os.mkdir(run_dir)
os.chdir(run_dir)
# generate regions:
region = []
rmax = wmax
while rmax > wmin:
	region.append((rmax - rmax * vc, rmax))
	rmax = rmax - rmax * vc2

# check last region isn't too small
last_max = region[-1][1]
region[-1] = wmin, last_max
# if it is too small, then join last and second last regions
if (last_max - wmin)/last_max < vc: 
	del region[-1]
	region[-1] = wmin, region[-1][1]
# first region (no need for including earlier fits)

wamin,wamax = region[0]

buffer = '%%%% ../%s   1   %10.5f  %10.5f vsig=%s\n\n' % (
	specfile, wamin, wamax, vsig)
open(inputfile,'w').write(buffer)
open('temp_vpfit_input_file.txt','w').write(
	'F ad %s il\n\n%s\nn\nn\n' % (chi_thresh,inputfile))
os.system('vpfit < temp_vpfit_input_file.txt')

badfit = []

# Choose most recent filename containing 'f26' as the output.
os.system('ls -lt > temporary_lsminuslt_list')
l = open('temporary_lsminuslt_list').readlines()
filename = [ line.split()[-1] for line in l ] 

for f in filename:
	if 'f26' in f:
		f26name = f
		break
# read in f26 file, and remove any lines which fell outside the
# fitting region, as they're likely to be wrong, and will be better
# fitted by other regions.  Also note if the chi-square threshold value
# was not reached
rows = open(f26name).readlines()
chisq = float(rows[1].split(None,4)[3])  # 4th item on second row
if chisq > chi_thresh:  badfit.append(f26name)
hd = rows[:2]
ions = rows[2:]
wav = [(float(ion[9:].split(None,1)[0]) + 1) * lya for ion in ions]
new_ions = [ions[i] for i in xrange(len(ions)) if wamin < wav[i] < wamax]
if len(new_ions) < len(old_ions):
	newrows =  hd + new_ions
	open(f26name,'w').writelines(newrows)
	if f26name not in badfit:  badfit.append(f26name)

# Now run vpfit on each region
for wamin,wamax in region[1:]:
	# Make an input (f26 format?) file
	output = ' %%%% ../%s   1   %10.5f  %10.5f vsig=%s\n\n' \
			 % (specfile, wamin, wamax, vsig)
	open(inputfile,'w').write(output)

	# Make the text file to give to vpfit (vpfit < infile).  Needs the
	# name of the input file and the F ad 2.0 il inc 'prevf.26 blah'
	# line.
	output = ["F ad %s il inc '%s 26 0 0 0 0 0 0 0 i'\n" %
			  (chi_thresh,f26name),"\n","%s\n" % inputfile,"n\n","n\n"]
	open('temp_vpfit_input_file.txt','w').writelines(output)
	os.system('vpfit < temp_vpfit_input_file.txt')

	# After vpfit finishes, read in the new f26 output file
	os.system('ls -lt > temporary_lsminuslt_list')
	l = open('temporary_lsminuslt_list').readlines()
	filename = [ line.split()[-1] for line in l ]
	for f in filename:
		if 'f26' in f:
			f26name = f
			break
	# read in f26 file, and remove any lines which fell outside the
	# fitting region, as they're likely to be wrong, and will be better
	# fitted by other regions.  Also note if the chi-square threshold value
	# was not reached
	rows = open(f26name).readlines()
	chisq = float(rows[1].split(None,4)[3])  # 4th item on second row
	if chisq > chi_thresh:  badfit.append(f26name + '\n')
	hd = rows[:2]
	ions = rows[2:]
	wav = [(float(ion[9:].split(None,1)[0]) + 1) * lya for ion in ions]
	new_ions = [ions[i] for i in xrange(len(ions)) if wamin < wav[i] < wamax]
	if len(new_ions) < len(old_ions):
		newrows =  hd + new_ions
		open(f26name,'w').writelines(newrows)
		if f26name not in badfit:  badfit.append(f26name + '\n')

# finished!  Search through f26 files and look for possibly dodgy
# fits, write these to a log file.
open('bad_fits.txt','w').writelines(badfit)
os.system('cat f26.* > all_f26')
os.chdir('..')
